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ABSTRACT 

Hierarchical Cold Dark Matter (CDM) models predict that Milky Way sized halos contain several 
hundred dense low-mass dark matter satellites (the substructure) , an order of magnitude more than the 
number of observed satellites in the Local Group. If the CDM paradigm is correct, this prediction implies 
that the Milky Way and Andromeda are filled with numerous dark halos. To understand why these halos 
failed to form stars and become galaxies, we need to understand their history. We analyze the dynamical 
evolution of the substructure halos in a high-resolution cosmological simulation of Milky Way sized halos 
in the ACDM cosmology. We find that about 10% of the substructure halos with the present masses 
< 10 8 — 10 9 M Q (circular velocities V m < 30 km/s) had considerably larger masses and circular velocities 
when they formed at redshifts z > 2. After the initial period of mass accretion in isolation, these objects 
experience dramatic mass loss due to tidal stripping. Our analysis shows that strong tidal interaction 
is often caused by actively merging massive neighboring halos, even before the satellites are accreted by 
their host halo. These results can explain how the smallest dwarf spheroidal galaxies of the Local Group 
were able to build up a sizable stellar mass in their seemingly shallow potential wells. We propose a 
new model in which all of the luminous dwarf spheroidals in the Local Group are descendants of the 
relatively massive (> 10 9 M Q ) high-redshift systems, in which the gas could cool efficiently by atomic 
line emission and which were not significantly affected by the extragalactic ultraviolet radiation. We 
present a simple galaxy formation model based on the trajectories extracted from the simulation, which 
accounts for the bursts of star formation after strong tidal shocks and the inefficiency of gas cooling in 
halos with virial temperatures T V1T < 10 K. Our model reproduces the abundance, spatial distribution, 
and morphological segregation of the observed Galactic satellites. The results are insensitive to the 
redshift of reionization. 

Subject headings: cosmology: theory-galaxies: formation-galaxies: dwarf-galaxies: halos- halos: 
evolution- methods: numerical 



1. INTRODUCTION 

Semi-analytic models of galaxy formation (Kauffmann 
et al. 1993; Bullock et al. 2000; Somerville 2002; Benson 
et al. 2002) and numerical simulations (Klypin et al. 1999b; 
Moore et al. 1999a) have convincingly showed that the ex- 
pected number of dark matter clumps around the galactic 
Milky Way (MW) sized halos exceeds the observed num- 
ber of satellites by an order of magnitude. The discrep- 
ancy may indicate that the amplitude of the small-scale 
primordial density fluctuations is considerably lower than 
expected in the Cold Dark Matter (CDM) scenarios (e.g., 
Kamionkowski & Liddle 2000; Zentner & Bullock 2003) 
or that dark matter is self-interacting (Spergel & Stein- 
hardt 2000). An alternative "astrophysical" interpreta- 
tion is that the mismatch indicates that galaxy formation 
in dwarf halos is inefficient. 

Several plausible physical processes may suppress gas 
accretion and star formation in dwarf dark matter (DM) 
halos. The cosmological UV background, which reion- 
ized the Universe at z > 6, heats the intergalactic gas 
and establishes a characteristic time-dependent minimum 
mass for halos that can accrete gas (e.g., Efstathiou 1992; 
Thoul & Weinberg 1996; Quinn et al. 1996; Navarro & 
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Steinmetz 1997; Gnedin & Hui 1998; Kitayama & Ikeuchi 
2000; Gnedin 2000; Dijkstra et al. 2003). The gas in the 
low-mass halos may be photoevaporated after reionization 
(Barkana & Loeb 1999; Shaviv & Dekel 2003; Shapiro 
et al. 2003). In particular, Shaviv & Dekel (2003) re- 
cently argued that halos with circular velocities of up to 
~ 30 kms -1 can be photo-evaporated by the UV back- 
ground. At the same time, the ionizing radiation may 
quickly dissociate molecular hydrogen, the only efficient 
coolant for low-metallicity gas in such halos, and prevent 
star formation before the gas is completely removed (e.g., 
Haiman et al. 1997). 

The combined effect of these processes is likely to leave 
all DM halos with masses < few x 10 9 M Q dark. This is 
consistent with current observational constraints which in- 
dicate that halos with M < 10 10 M Q are virtually devoid 
of galaxies (van den Bosch et al. 2003). It is thus remark- 
able that the dynamical masses of some of the Local Group 
dwarfs are only ~ 10 7 M (Mateo 1998). How could such 
galaxies form stars despite the suppressing processes listed 
above? 

One possibility is that they manage to accrete a certain 
amount of gas before the Universe is reionized (Bullock 
et al. 2000) with the implicit assumption that this gas 
can be subsequently converted to stars. However, it is 
likely that gas cooling and star formation in such small 
systems is inefficient. For example, cosmological simula- 
tions with self-consistent treatment of H2 chemistry and 
radiative transfer indicate that star formation is strongly 
suppressed in halos with masses M < 5 x 10 8 M Q at all red- 
shifts, even before reionization (Chiu et al. 2001). In addi- 
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tion, the galaxies may not be able to form sufficiently early 
to accrete the gas in the first place, if the power spectrum 
normalization is low or the Universe was reionized early, as 
indicated by the first-year WMAP results (Spergel et al. 
2003; Kogut et al. 2003). An alternative proposal was re- 
cently suggested by Stoehr et al. (2002, 2003) and corrab- 
orated by Hayashi et al. (2003), who argued that the host 
halos of the low-luminosity dwarf spheroidal galaxies may 
be considerably more massive than previously thought. In 
this case, the large halo mass could allow an object to 
resist the suppressing effects of UV background. 

In this paper we study the dynamical evolution of dwarf 
satellite halos around the Milky Way-sized hosts in self- 
consistent cosmological simulations. We show that the 
evolution of such objects is complex and often involves 
dramatic tidal stripping, interactions with other satellites, 
mass loss, and changes of internal structure. Most im- 
portantly, we find that some of the satellites that have 
small masses and circular velocities at the present, once 
were considerably more massive and could have plausibly 
formed stars in the past. We argue that the evolution of 
these objects may explain how the smallest dwarfs in the 
Local Group managed to form their stellar populations. 

The paper is organized as follows. In § 2 we describe 
the details of the numerical simulation used in our anal- 
ysis. In § 3 and § 4 we discuss the algorithm used to 
identify halos and the method used to construct their evo- 
lutionary tracks. In § 5 we present the main results on the 
dynamical evolution, abundance, and radial distribution 
of the dwarf dark matter halos. In § 6 we present a model 
for star formation in these systems and compare results 
to the observed abundance and spatial distribution of the 
Local Group dwarfs. We discuss the implications of our 
results and compare our model to the previous studies in 
§§7 and 8. Finally, in § 9 we summarize our findings and 
conclusions. 

2. SIMULATION 

We used the Adaptive Refinement Tree TV-body code 
(ART, Kravtsov et al. 1997; Kravtsov 1999) to follow the 
evolution of three Milky Way-sized halos in the concor- 
dance ACDM cosmology: (n m ,Q A ,h,a 8 ) = (0.3,0.7,0.7,0.9) 
The simulation starts with a uniform 256 3 grid covering 
the entire computational box. This grid defines the low- 
est (zeroth) level of resolution. Higher force resolution is 
achieved in the regions corresponding to collapsing struc- 
tures by recursive refining of all such regions using an 
adaptive refinement algorithm. Each cell can be refined 
or de-refined individually. The cells are refined if the par- 
ticle mass contained within them exceeds a certain speci- 
fied value. The grid is thus refined to follow the collapsing 
objects in a quasi-lagrangian fashion. 

The galactic halos were simulated in the comoving box of 
25/i _1 Mpc; they were selected to reside in a well-defined 
filament at z = 0. Two halos are neighbors, located at 
425/i _1 kpc (i.e., 610 kpc <~ 2i? v j r ) from each other. 
The configuration of this pair thus resembles that of the 
Local Group. The third halo is isolated and is located 
<~ 2 Mpc away from the pair. 

Multiple mass resolution technique was used to set up 
initial conditions. Namely, a lagrangian region correspond- 
ing to a sphere of radius equal to two virial radii around 



Table 1. Properties of the host halos 



Halo M vir R vir V m Environment 

(/i^M©) {h^ 1 kpc) (kms- 1 ) 



Gi 


1.66 x 10 12 


298 


213 


isolated 


G 2 


1.24 x 10 12 


278 


199 


pair 


G 3 


1.19 x 10 12 


281 


183 


pair 



Note - i? v i r is the virial radius corresponding to the aver- 
age density of 180 times the mean density of the universe in 
h^ 1 kpc; M v i r = M(< R v i T ) in ft _1 M (both radius and mass 
are given for z = 0); V m is the maximum circular velocity. 

each halo was re-sampled with the highest resolution par- 
ticles of mass 77j p = 1.2 x lO 6 /i _1 M , corresponding to 
1024 3 particles in the box, at the initial redshift of the 
simulation (zj = 50). The high mass resolution region 
was surrounded by layers of particles of increasing mass 
with a total of 5 particle species. Only regions containing 
highest resolution particles were adaptively refined and the 
threshold for refinement was set to correspond to the mass 
of the four highest resolution particles. The maximum of 
ten refinement levels was reached in the simulations cor- 
responding to the peak formal spatial resolution of 150 
comoving parsec. Each host halo is resolved with ~ 10 6 
particles within its virial radius at z = 0. 

From this point, we will refer to the isolated halo as Gi 
and the halos in pair as G 2 and G 3 . These halos are called 
Bi, Ci, and Di, respectively, in Klypin et al. (2001) and we 
refer the reader to this paper for further details. The main 
properties of these three host halos, the virial mass, radius, 
and maximum circular velocity, are given in Table 1. We 
choose to define the virial radius (and the corresponding 
virial mass) as the radius encompassing the density of 180 
times the mean density of the universe. For the commonly 
used overdensity of 340, the virial radii and masses for Gi, 
G 2 , and G3 arc i? 3 4 = 231, 2 1 2, and 213h~ 1 kpc and 
M 340 = 1.45 x 10 12 , 1.13 x 10 12 , and 1.14 x lO 12 ^ 1 M , 
respectively. The masses are in the range of possible MW 
halo masses (Klypin et al. 2002). 

Figure 1 shows the mass aggregation history of the three 
host halos. They have similar masses at the present but 
rather different evolutionary histories. In all cases, there is 
a period of very rapid mass assembly at z > 2 — 3 followed 
by a relatively quiescent accretion at z < 1.5, the behav- 
ior typical of hierarchically forming halos (Wechsler et al. 
2002). Host Gi undergoes a spectacular multiple major 
merger at z ss 2, which results in a dramatic mass increase 
on a dynamical time scale. Halos G2 and G3 increase their 
mass in a series of somewhat less spectacular major merg- 
ers which could be seen as mass jumps at 5 < z < 1. All 
three systems accrete little mass and experience no major 
mergers at z < 1 (or lookback time of w 8 Gyr) and thus 
could host a disk galaxy. Note, however, that halos Gi 
and G3 experience minor mergers during this period. 

3. halo identification 

In this study we use a variant of the Bound Density 
Maxima (BDM, Klypin et al. 1999a) halo finding algo- 
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Fig. 1. — Mass aggregation histories for the three MW-size host 
halos analyzed in this study. 



rithm to identify halos both within (subhalos) and outside 
the host halos. Throughout this paper we will use the 
terms subhalo, substructure, and satellite interchangeably 
to indicate the distinct gravitationally self-bound halos lo- 
cated within the virial radius of a larger halo, which we 
call the host. The division is illustrated in Figure 2. 

The BDM algorithm first finds positions of local maxima 
in the density field smoothed on a certain scale. Starting 
with the highest overdensity particle, we surround each 
potential density maximum by a sphere of radius r^ n( \ = 
lOh^ 1 kpc and exclude all particles within this sphere from 
further search. The search radius is defined by the size of 
smallest systems we aim to identify. We verified that the 
results do not change if this radius is decreased by a factor 
of up to four. After all potential halo centers are iden- 
tified, we analyze the density distribution and velocities 
of surrounding particles to test whether the center corre- 
sponds to a gravitationally bound clump. Specifically, we 
construct the density, circular velocity, and velocity dis- 
persion profiles around each center and iteratively remove 
unbound particles (see Klypin et al. 1999a, for details). 
We then construct final profiles using only bound parti- 
cles and use them to calculate such halo properties as the 
maximum circular velocity V m , mass M , etc. 

The virial radius is meaningless for the subhalos within a 
larger host as their outer layers are tidally stripped and the 
extent of the halo is truncated. The definition of the outer 
boundary of a subhalo and its mass are thus somewhat am- 
biguous. We adopt the truncation radius, r t , at which the 
logarithmic slope of the density profile constructed from 
the bound particles becomes larger than —0.5 as we do 
not expect the density profile of the CDM halos to be 
flatter than this slope. Empirically, this definition roughly 
corresponds to the radius at which the density of the grav- 
itationally bound particles is equal to the background host 
halo density, albeit with a large scatter. For some halos 




Fig. 2. — Distribution of dark matter particles (points) and dark 
matter halos (circles) identified by our halo finding algorithm cen- 
tered on the isolated galactic halo at z = 0. The radius of the 
largest circle indicates the actual virial radius, R vlr , of the host halo 
(i? v ir = 298/i _1 kpc); the radii of the other halos are the minimum 
between truncation radius rt and i? v ir- The particles are colored 
on a gray-scale logarithmic stretch according to their local density. 
The stretch is chosen to highlight the cores of the halos for clarity. 



r t is larger than their virial radius. In this case, we set 
r t — i? v ir- Throughout this paper, we will denote the min- 
imum of the virial mass and mass within r t , simply as 
M. For each halo we also construct the circular velocity 
profile V c (r) — GM (< r)/r and compute the maximum 
circular velocity profile V m . 

Figure 2 shows the particle distribution in the halo Gi 
at z = along with the halos (circles) identified by the 
halo finder. The particles are color-coded on a gray scale 
according to the logarithm of their density to enhance vis- 
ibility of substructure clumps. The radius of the largest 
circle indicates the actual virial radius, i? v ir, of the host 
halo (i? V ir = 298/i~ 1 kpc); the radii of the other halos are 
the minimum of the truncation radius r t and i? V ir- The 
figure demonstrates that the algorithm is efficient in iden- 
tifying the substructure down to small masses. 

4. CONSTRUCTING TRAJECTORIES 

The halo finder described above was run at the 96 saved 
epochs between z = 10 and z = with a typical spacing of 
~ 1 — 2 x 10 8 yr between outputs. For each epoch, the halo 
finder produced a halo catalog with positions, velocities, 
radii i\ = min(rt, r v i r ), masses m(< r^), maximum of the 
circular velocity profile V m and the radius at which the 
maximum occurs r max . In addition, for each halo we save 
indices of all gravitationally-bound DM particles located 
within Th. 

This information is used to identify the progenitors of 
halos at successive epochs. Specifically, for a current epoch 
Zi, starting at z = 0, we search progenitors for each halo 
at several previous epochs Zi_j as follows. First, we select 
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a given fraction, /bound, of the most bound particles of the 
halos at the epochs of consideration. We then compare 
the fraction of these particles that is common between all 
pairs of halos at successive epochs and assume that the 
halo with the highest common fraction is the progenitor. 
The trajectories used in this study were constructed using 
/bound = 0.25. As the halo catalogs may miss some halos, 
especially near the completeness limit of the simulation, 
if the progenitor is not found at the previous epoch we 
need to search at the earlier epoch, etc. In particular, if 
a halo is located within the search radius r soarcn of some 
larger system it will not be identified by the halo finder. In 
constructing the trajectories we search for progenitors of a 
halo at z\ at epochs up to z;_4. In the dominant majority 
of cases the progenitors are found at the previous epoch 
z;_i. We experimented with other algorithms for progen- 
itor identification and found the adopted prescription to 
be the most reliable and efficient. 

5. RESULTS 

5.1. Tidal stripping and dynamical evolution of satellite 

halos 

Figure 3 shows three examples of the evolution of satel- 
lite halos. In the middle row of panels we plot the tidal 
force experienced by each object. The force was calculated 
both directly from the gravitational potential field com- 
puted in the simulation and analytically from the neigh- 
bor halo catalogs, as described in Appendix A. The Figure 
shows the trace of the tidal tensor, F t ;d = J2 a Faa, which 
is a good measure of the overall tidal field. 

In all cases the tidal force experienced by the satellite 
coincides approximately with the time when the object is 
closest to the host, as expected. The Figure shows, for 
example, that at later epochs the tidal force calculated 
directly using the potential from the simulation can be 
well approximated by the analytical force from the host 
halo (see eq. A4). However, at earlier epochs (e.g., the 
highest peak in the left column) the force from the host 
underestimates the total tidal force. Thus, the overall tidal 
stripping is produced not only by the host halo but also by 
the massive neighbor halos, even before the host is formed 
(Gnedin 2003b). 

The tidal heating by multiple halos is similar to "galaxy 
harassment" in clusters of galaxies (Moore et al. 1996, 
1999b), except that it may occur when the halo is still 
isolated. As Figure 3 shows, the true force computed from 
the potential can be recovered if the analytical contribu- 
tions of all neighboring halos are included. Their contri- 
bution is particularly important during major mergers of 
the host, when the centers of two or more massive halos 
are located in the close vicinity of each other and the satel- 
lite halos. Our analytic estimate describes the strong tidal 
peaks remarkably well, but becomes inaccurate for low (a 
few Gyr" 2 ) values of F t id- 

The amount of energy imparted to the halo depends on 
the square of the tidal force (eq. A6). Thus, by far the 
strongest tidal heating experienced by an object is during 
the highest tidal peaks. For the object in the left column 
of Figure 3, for example, most of the stripping and dis- 
ruption is due to the tidal peak at t ~ 4 Gyr (z ~ 1.5). 
At this epoch, the host halo is not yet fully assembled and 
is undergoing a major merger with three other massive 



halos. It is at this epoch, however, that the satellite ex- 
periences the most dramatic tidal mass loss. Subsequent 
tidal peaks result only in a relatively mild stripping. The 
object in the middle panel also suffers a dramatic mass loss 
at t ~ 4 — 5 Gyr. In this case, however, the efficient tidal 
stripping continues due to the later pericentric passages 
and associated peaks in the tidal force. Finally, the third 
satellite shown in the figure experiences only a relatively 
mild tidal stripping. This satellite orbits in the outer re- 
gions of the host and never reaches the central w 60 kpc. 

Note that the pericenter of the third satellite is larger at 
the late epochs compared to the pericenter at t w 4 Gyr. 
This is contrary to a naive expectation that the pericenter 
should stay constant or decrease with time if dynamical 
friction is efficient. The real situation is clearly more com- 
plicated. The satellite can lose as well as gain the orbital 
energy. The latter can occur via a three-body interaction. 
Indeed, in examining individual trajectories we found cases 
where a satellite gains orbital energy via the "slingshot" 
acceleration — a classic three-body interaction. 

Figure 3 demonstrates that some satellites with small 
maximum circular velocity and mass at z = were sub- 
stantially more massive during the early stages of their 
evolution. The mass of the object in the middle column 
is i=a 10 10 M Q and its circular velocity is > 40 kms -1 at 
t = 4 Gyr. At the present epoch, they are only 2 x 10 8 M 
and 18 kms -1 , respectively. In the extreme cases we find 
changes of mass and V mayi by a factor of 200 and 8, respec- 
tively (see Fig. 6). 

At the same time, the object in the right column of Fig- 
ure 3 has a considerably larger pericenter and experiences 
weaker tidal force by more than an order of magnitude. 
Consequently, its mass and circular velocity change little 
during the evolution. What is the relative frequency of 
such cases compared to the cases of dramatic mass loss? 
We address this question in the next section. 

5.2. Internal structure evolution 

Strong tidal forces experienced by orbiting halos lead to 
a substantial mass loss, preferentially at the outer radii. 
The changes in the inner regions are more subtle and occur 
at a slower rate but can nevertheless be significant (e.g., 
Klypin ct al. 1999a; Hayashi et al. 2003; Stoehr et al. 2003; 
Kazantzidis et al. 2003). Figure 4 shows the maximum 
values of M and V max reached by a satellite during its 
evolution vs. their present values. Most of the surviving 
satellites experience only mild evolution, less than a factor 
of two in V m . Yet, there is a fair number of cases in which 
the evolution is significant. The average changes in M or 
V m do not seem to depend on the halo mass. 

Figure 5 shows the ratio of the mass at z — to the 
maximum mass achieved by each satellite during its evo- 
lution vs. the ratio of the maximum circular velocities at 
these two epochs. The figure shows a strong correlation 
between the two ratios: 

j^max = ( ymax ) ' ~ ^ - 4, (1) 

where and M° are the values at the present and V™ ax 
is the maximum circular velocity at the epoch when the 
halo reached the maximum mass, M max . This correlation 
shows that the internal structure of satellites re-adjusts 
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t (Gyr) t (Gyr) t (Gyr) 

Fig. 3. — Three examples of the evolution of satellites of a MW-size host (different columns). Top panels: the proper distance between the 
satellite and the center of the host halo as a function of time. Middle panels: the tidal force experienced by the object, calculated directly 
from gravitational potential, is shown by the solid line. Dotted line shows the equivalent tidal force from the host halo with the host density 
profile approximated by an NFW model. Dashed line shows the contributions of all neighboring halos, including the host, with the density 
profiles of halos approximated by an NFW model with r max and Vm as measured in the simulation. See Appendix A for details on the tidal 
force calculation. Bottom panels: maximum circular velocity V m (solid lines) and bound mass m(< rt) (dashed lines) as a function of time. 
The three objects show different types of evolution: dramatic early stripping with a relatively quiescent evolution afterward (left), continuous 
dramatic tidal stripping (middle), weak stripping and quiescent evolution (right). 
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Fig. 4. — Top panel: the maximum circular velocity along the 
trajectory, V^ 3 -", vs. the present-day value VC2 . Bottom panel: 
maximum mass along the trajectory vs. the present-day mass of a 
halo. The satellite halos shown are located within the virial radius 
of their respective host halo. The figure shows that many halos 
experience a dramatic decrease in their mass and circular velocity. 
The solid circles show the halos that host luminous galaxies in our 
model (see § 6). 



as they lose mass due to tidal stripping. Note that the 
decrease of V m indicates the decrease in density within the 
inner radius of ss 2.16r s , where r s is the NFW scale radius. 
The adjustment is such that the virial correlation 



MocV", a«3 



(2) 



is approximately maintained at all times. 

This can be seen in Figure 6, which shows the tracks 
of individual satellites in the M — V m plane. The satel- 
lites shown were selected from all three galactic hosts. We 
selected objects with large changes in mass to maximize 
the dynamic range. The figure shows that both during the 
periods of mass growth while evolving in isolation and the 
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Fig. 5. — The ratio of the 2 = mass to the maximum mass 
achieved by each satellite during its evolution vs. the ratio of the 
maximum circular velocities at these two epochs. The dots represent 
the ratios for individual halos. Solid circles show the average for the 
equally space logarithmic mass bins and solid line shows the power 
law weighted least square fit to these points. Open circles and the 
dashed line show the same for the binning in V^n ax - The histograms 
in the top and right panels show the fraction of the halos with a 
given mass and circular velocity ratio in logarithmic bins of size 0.2 
and 0.05 for the mass and velocity ratios, respectively. 



periods of mass decrease due to tidal stripping, halos ap- 
proximately move up and down the power-law dependence 
of eq. (2). For instance, the track shown by in the top left 
panel starts at M « 3 x 10 8 /i _1 M© and V m w 30 kins" 1 
at z = 10 and by the redshift z = 2 reaches the mass of 
2 x 10 10 /i _1 M©. In the ten billion years between z = 2 
and z = 0, the halo loses 99.5% of the mass and its V m de- 
creases by a factor of eight. Yet, during the entire course 
of evolution the halo moves roughly along the M oc V^ 3 
line. 

This result is in agreement with Hayashi et al. (2003), 
who found a correlation similar to that of eq. (1) using con- 
trolled TV-body experiments to study the tidal stripping 
and internal structure of the NFW halos (see their Fig. 
12). They note that the density at all radii changes in re- 
sponse to tidal shocking. The density decrease is greatest 
at large radii, so that the overall profile steepens while the 
normalization drops. The adjustment of the density profile 
leads to the decrease of both V m and r max . We also find 
that r max of satellites in our simulations decrease system- 
atically as they lose mass. A similar evolution of V max as a 
function of mass loss is found in very high-resolution con- 
trolled simulations of Kazantzidis et al. (2003), which fol- 
lowed tidal stripping of an NFW satellite resolved with 10 7 
particles (Stelios Kazantzidis, private communication). 
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Fig. 6. — Trajectories (solid lines) of halos on the mass — max- 
imum circular velocity plane for four alos found within the virial 
radius of all three galactic hosts at z = 0. We selected objects 
with large changes in mass to maximize the dynamic range and il- 
lustrate the effect. The solid circles mark the end (z = epoch) 
of each trajectory. The figure shows that both during the periods 
of mass growth while evolving in isolation and the periods of mass 
decrease due to tidal stripping, halos approximately move up and 
down the virial dependence M oc V£ with a ~ 3 — 5; the dotted line 
corresponds to a = 3.3. 



5.3. Evolution of halos in the M — V m plane 

In the previous section we showed that individual halos 
maintain M oc V£ relation during their evolution. This 
explains why the same correlation between the mass and 
V m holds for both subhalos and isolated halos (Avila-Reese 
et al. 1999; Bullock et al. 2001). We also find that the 
mass-circular velocity relations for the subhalos at z = 
and for their progenitors at the epoch when the maximum 
mass was reached have similar amplitudes and slopes (w 
3.3). 

In this section we will consider the mechanism behind 
such behavior in more detail. We can interpret the ob- 
served evolution of mass and V" m for a given subhalo by 
dividing it in the following two stages: (1) mass growth 
while evolving in isolation, and (2) mass decrease due to 
tidal stripping after the halo is accreted by its host. The 
transition between the two stages typically occurs at z ~ 2 
for the mass range of subhalos and hosts considered here. 

We fit the slope a for the trajectories in the M — V m 
plane for all satellites separately in the two regimes. We fit 
only tracks of halos with V m > 15 km s" 1 at z — and with 
at least 10 rcdshift outputs. In calculating the average 
slope, an, and the dispersion of the sample, a aj we weigh a 
of each halo by the error of the fit. The isolated halos have 
the average slope a = 4.7, with the dispersion a a = 1.2. 



The truncated halos have an = 2.9 and o~ a = 1.2. Thus the 
slopes in the two regimes seem to be somewhat different. 

These different slopes can be linked to the different av- 
erage densities of the dwarf halos in the two regimes. In 
the mass-growth stage, when the average density of the 
Universe is p(z) = p (l + z) 3 , the mass and velocity are 
given by the virial scaling relation M v ; r oc V^ 3 ir (l + z)~ 3 / 2 . 
Also, Bullock et al. (2001) showed that, as long as the 
NFW model is an adequate description of the halo pro- 
file, V" m and Kir are related through the concentration 
parameter approximately as V m /V v i T oc cV 4 . The me- 
dian concentration itself varies with the mass and red- 
shift as c v i r oc (1 + z)^ 1 M~® A3 . Since all the variables 
scale as some power of (1 + z), it is natural to approxi- 
mate the evolution of the mass and maximum velocity as 
M v ir oc (1 + z)~ q and V m oc (1 + z)~ p , with the above re- 
lations leading to ||g — 3p + |. The slope an = q/p = 4.7 
is achieved for q — 2.8 and p = 0.6, although the scatter 
in the value of a implies a corresponding scatter in the 
exponents q and p. 

Note that the slope of the M v ; r — V m relation is steeper 
than the virial a w 3 because the virial parameters of 
isolated halos depend on the mean density of the Universe 
and that density is changing with redshift. The same zero- 
point of the relation can be maintained only if both M v ; r 
and V m are changing with time in a certain way, specified 
above. 

During the second stage of evolution, the subhalo expe- 
riences tidal forces from the host and other halos, and its 
mass and extent are tidally truncated. The average density 
p t within the truncation radius i?t is approximately con- 
stant along the orbit and is proportional to the background 
density of the host halo at the pericenter of subhalo's orbit. 
The truncation radius scales with the truncated mass, M t , 

as R t oc (M t /p t ) 1/3 , so that M t oc V t 3 p^ 1/2 . The velocity 
V t oc (Mt/Rt) 1 / 2 is a good estimate for the peak velocity 
when the subhalo is severely truncated. If the background 
density is constant along the satellite trajectory, we ob- 
tain the following relation: M t oc V^, in agreement with 
the average fit in this regime. Figure 6 shows that in this 
regime the power-law relation between M t and V m has a 
significant dispersion, which is due to the variation of pt 
along the trajectory. Nevertheless, as long as the distance 
of closest approach to the host halo remains the same, the 
average relation is well maintained. 

Thus, we expect the mass- velocity relation to be con- 
strained by the two limiting slopes 3 and 5. The actual 
mass accretion and mass loss history may vary from halo to 
halo but the same M — V m relation is maintained through- 
out the evolution, with a transition from the initial slope 
a = 4.7 for the isolated halos to the later slope a = 3 for 
the tidally truncated halos. 

5.4. Abundance and radial distribution of galactic 
satellites 

Figure 7 shows the cumulative velocity functions (CVFs), 
the number of satellites with maximum circular velocity 4 

4 Note that uncertainty in the velocity anisotropy affects the conver- 
sion of the line-of-sight rms velocity of dSph galaxies to Vm. In the 
plot we assume an isotropic velocity distribution. Our re-analysis of 
numerical simulations of Gnedin (2003a) shows that tidal truncation 
and heating of galaxies leads to the preferential removal of radial 
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Fig. 7. — The cumulative velocity function of the dark matter 
satellites in the three galactic halos (solid lines compared to the 
average cumulative velocity function of dwarf galaxies around the 
Milky Way and Andromeda galaxies (stars). For the objects in 
simulations V c i rc is the maximum circular velocity, while for the Lo- 
cal Group galaxies it is either the circular velocity measured from 
rotation curve or from the line-of-sight velocity dispersion assum- 
ing isotropic velocities. Both observed and simulated objects are 
selected within the radius of 200h~ 1 kpc from the center of their 
host. The dashed lines show the velocity function for the luminous 
satellites in our model described in § 6. The minimum stellar mass 
of the luminous satellites for the three hosts ranges from rs 10 5 Mq 
to f» 10 6 Mq, comparable to the observed range. 



larger than a given value, for the objects located within 
200ft- 1 kpc of their host halo. The figure compares the 
CVFs for the DM satellites and observed satellites of the 
MW and Andromeda 5 and highlights the "missing satel- 
lite problem" (Kauffmann et al. 1993; Klypin et al. 1999b; 
Moore et al. 1999a): a large difference in the number of 
dwarf-size DM satellites in simulations and the observed 
number of dwarfs in the Local Group. 

Figure 8 shows the normalized cumulative radial distri- 
bution of the DM satellites compared to the radial distri- 
bution of satellites around the Milky Way within the same 
radius. The Local Group data is from the compilation of 
Grebel et al. (2003). The figure clearly shows that the spa- 
tial distribution of dwarf galaxies around the Milky Way 
is more compact than the distribution of the DM popu- 

orbits and the development of the tangentially-biased dispersion in 
the outer parts. A similar result has been found by Kazantzidis et al. 
(2004) and Moore et al. (2003). The solution of the Jeans equation 
for Vm is sensitive to the exact value of the anisotropy parameter 
(Zentner & Bullock 2003; Kazantzidis et al. 2003). 

5 We use the circular velocities compiled by Klypin et al. (1999b) 
with updated values of circular velocity for the Large and Small 
Magellanic Clouds of V m = 50 kms -1 and 60 kms" 1 , respectively 
(van der Marel et al. 2002) 



Fig. 8. — The fraction of satellites within a certain distance from 
the center of their host galaxy. The solid lines show distributions of 
the ACDM satellites in the three galactic halos, while the connected 
stars show the distribution of dwarf galaxies around the Milky Way. 
The figure shows that radial distribution of observed satellites is 
more compact than that of the overall population of dark matter 
satellites. The dashed lines show distributions for the luminous 
satellites in our model (§ 6). The population of luminous satellites 
is the same in this and previous figures. 



lation. The median distance of observed satellites within 
200ft,- 1 kpc is 60ft- 1 kpc and 85ft- 1 kpc for the MW and 
M31, respectively. For the DM satellites the correspond- 
ing median distances are 116ft- 1 kpc, 121ft- 1 kpc, and 
120ft- 1 kpc. Although the median for M31 satellites is 
smaller than that of the DM satellites, their radial distri- 
butions are formally consistent. However, the comparison 
with the M31 satellites is difficult at present because typ- 
ical distance errors are ~ 20 — 50 kpc (and > 70 kpc for 
some galaxies), comparable to the distance to the host. 

For the MW satellites the typical distance errors are an 
order of magnitude smaller and the comparison is consid- 
erably more meaningful. The Kolmogorov-Smirnov (KS) 
test gives probability of (6 — 8) x 10~ 4 that the MW satel- 
lites are drawn from the same radial distribution as the 
DM satellites. This has also been pointed out recently by 
Taylor et al. (2003), who compared the spatial distribution 
of the MW satellites to results of their semi-analytic model 
of galaxy formation. Thus, in addition to the vastly dif- 
ferent abundances of the observed and predicted satellites, 
there is a discrepancy in the radial distribution. Models 
that aim to reproduce the abundance of the LG satellites 
should therefore be able to reproduce the radial distribu- 
tion as well. 

6. A MODEL OF STAR FORMATION IN DWARF HALOS 
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6.1. Description of the Model 

In order to gain insight into which halos might become 
luminous and which might not, we implement the follow- 
ing simple model of star formation. We use the standard 
assumption that the gas within the halos with the virial 
temperature T v ; r > 10 4 K dissipates its energy via radia- 
tive cooling and forms a disk. We then apply the empirical 
Schmidt law to calculate the star formation rate in radial 
shells within the disk. The novel features of our model 
include: (i) use of mass accretion and stripping history 
of the dwarf halos extracted from simulation, (ii) effects 
of photoionizing extragalactic background using the filter- 
ing mass, (iii) effects of inefficient dissipation of the gas 
at Tvir ^ 10 4 K, and (iv) bursts of star formation due 
to strong tidal shocks. The details of the model are as 
follows. 

(i) Using the mass assembly history (MAH) of a given 
satellite halo directly from the simulation, instead of a 
semi-analytic approach, we are able to trace the major 
merger events as well as the quiescent accretion of ma- 
terial. The halo mass increases in both regimes, but the 
star formation rates are very different. The use of simula- 
tion MAHs allows us to determine the accretion epoch of 
a satellite and follow its mass loss due to tidal stripping. 

(ii) At each time output we calculate the accreted mass 
since the last time step, AM. We increase the total gas 
mass, M g , in the satellite by the amount of cold gas in 
a single halo with the mass AM at that epoch: AM S = 
f g (M, z)AM. The fraction f g takes into account the pho- 
toevaporation of baryons by extragalactic UV flux, using 
the filtering scale parametrization of Gnedin (2000) and 
taking the redshift of reionization to be z r = 7. Sec Ap- 
pendix B and equation (B3) for details. After the satellite 
enters the host halo, the accretion of new gas is halted and 
the disk scale length is fixed, although stars may continue 
to form from the remaining reservoir of cold gas. 

We distribute the gas on a spherically symmetric grid 
of 50 radial shells, according to the surface density of an 
exponential disk: £ 3 (r) = So exp (— r/r&). We use the ob- 
served Schmidt law of star formation to estimate the star 
formation rate: £* = 2.5x 10~ 4 (S g /M pc" 2 ) 1 - 4 M kpc~ 2 . 
Only the shells above the threshold S g > £ t h = 5 M Q pc~ 2 
form stars (Kennicutt 1998). 

(iii) The scale length of the disk is determined by its 
angular momentum. For a rotationally-supported disk it 
is approximately r<j = 2~ 1//2 A r v ; r . The value of the an- 
gular momentum parameter is drawn randomly from the 
probability distribution 
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with A = 0.045, o\ = 0.56, according to the latest mea- 
surement by Vitvitska et al. (2002). This is a key assump- 
tion of the semi-analytic models of galaxy formation. 

However, small halos at high redshift could only cool by 
atomic hydrogen to about 10 4 K. If their virial tempera- 
ture is only slightly above that equilibrium temperature, 
the gas would not be able to dissipate enough to reach 
a rotationally-supported state. Instead, its distribution 
would be more extended, which can have important impli- 
cation for the star formation with a density threshold £ t h- 
This effect is particularly important for dwarf halos. 



We model the effect of inefficient dissipation by adopt- 
ing the expansion factor that depends on the ratio of the 
virial temperature to the equilibrium temperature of 10 4 
K. The gas would reach a Boltzmann distribution with 
the density M/r 3 oc exp (— $/fcT), where <I> is the poten- 
tial energy Using the maximum circular velocity instead 
of the temperature and ignoring the slow variation of the 
potential, we can express the scale length of the gas as 
rd oc exp [c(V4/U m ) 2 ], where c is a normalization factor 
and V4 = 16.7 km s _1 is the virial velocity correspond- 
ing to T vir = 10 4 K. We find that c = 10 is a best fit to 
the abundance and radial distribution of the Local Group 
galaxies (see §6.2). This scaling also provides a good de- 
scription of the extent of the gas within halos in cosmolog- 
ical galaxy formation simulation described in Kravtsov & 
Gnedin (2003). Thus, we set the size of the gaseous disk 
at each time step to be 

r d = 2- 1 / 2 Xr vir xe 10 ^/ y ^ 2 . (4) 
Of course, is not allowed to exceed the tidal radius of 
the halo, r t . The gas in large halos with V m 3> V4 can cool 
efficiently and reach rotational support, but for small halos 
with V m > V4 the extended distribution reduces the central 
concentration of the gas and hinders star formation. 

(iv) Strong tidal forces, such as in the interacting or 
merging galaxies, may lead to a burst of star formation 
throughout the dwarf galaxy. The association of star- 
bursts with strong peaks of the tidal force is motivated 
by theoretical models (Mayer et al. 2001a) and, to a cer- 
tain extent, by observations (Zaritsky & Harris 2003). The 
latter suggest that the tidally-triggered star formation in 
the SMC can be accurately modeled as an instantaneous 
burst of star formation. Zaritsky & Harris (2003) find the 
best fit to their data when the star formation rate (SFR) 
varies as r~ 4 - 6 with the distance to the Galaxy. The tidal 
interaction parameter, J t id (see eq. [A7]), that reflects the 
integrated effect of a single tidal shock, is the most natural 
candidate for the parametrization of the tidally-triggered 
SFR. Ignoring the adiabatic correction, it varies with the 
distance to the perturber approximately as 7 t id oc r~ 4 (but 
see the discussion in § 6.2). 

We allow for the starburst mode of star formation, when 
the tidal interaction parameter exceeds a threshold value. 
After experimenting with different thresholds, we find that 
Im,th = 4x 10 3 Gyr -2 provides the best simultaneous fit to 
the velocity function and spatial distribution of the satel- 
lites. In all radial shells, a fraction/* = itid/4 x 10 4 Gyr~ 2 
(with a maximum of /* = 0.5) of the available gas is con- 
verted into stars instantaneously. The normalization of /* 
is somewhat arbitrary and can be adjusted to fit the stel- 
lar masses of the satellites. Since the starburst changes 
drastically the distribution of gas in the galaxy, new in- 
falling gas may have a very different angular momentum. 
Therefore, after each starburst we recalculate the value of 
A according to eq. (3). 

The external tidal force determines the truncation ra- 
dius R t of the satellite, outside which all stars and gas are 
lost. In a static gravitational field, the radius of the Roche 
lobe is set by the condition that the average density of mat- 
ter in the satellite equals twice the local ambient density 
(for the isothermal sphere potential). In a dynamic sit- 
uation of the satellite on an eccentric orbit experiencing 
tidal shocks, the truncation depends on the time-varying 
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tidal force. However, using iV-body simulations of the dy- 
namical evolution of galaxies in clusters, Gnedin (2003a) 
showed that the truncation radius can be accurately de- 
scribed by the same condition, p a - v {Rt) — 2 p t id, where the 
effective tidal density ptid is related to the trace of the 
tidal tensor via 

PtM - 1.8 x 10- 5 ((Tpa) Af QP c- 3 . (5) 

The truncation occurs near the maximum of the tidal force 
along the orbit, usually at the perigalactic distance. 

The knowledge of the external tidal force also allows us 
to estimate the tidal heating of stars in the satellite. After 
each tidal shock, typically once per orbit, the velocity dis- 
persion of stars in each radial shell increases by the amount 
(Gnedin 2003b) 



Acr 2 (r) = 0.32 
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The mass- weighted dispersion a may serve as an indicator 
of the morphological type of the satellite. In §6.2, we adopt 
the ratio of the rotation velocity to the velocity dispersion, 
v ro t/a, as a possible criterion. In practice, we compute 
Urot as the circular velocity of the NFW halo at the radius 
enclosing all bound stars. 

6.2. Results 

We show the predictions of our model for the cumula- 
tive velocity function and radial distribution of luminous 
satellites by dashed lines in Figures 7 and 8. The model 
reproduces fairly well both observational statistics. The 
abundance of luminous satellites and the shape of the ve- 
locity function are in reasonable agreement with observa- 
tions. The stellar masses of the satellites are in the range 
10 5 < M* < 10 10 Mq, similar to the observed range of stel- 
lar masses of the LG galaxies (e.g., Dekel & Woo 2003). 
As can be seen in Figure 4, all of the luminous systems in 
our model, including those that have masses and circular 
velocities of the smallest dwarfs, were relatively massive 
(M > 10 9 M and V m > 30 kms" 1 ) at some point in 
their evolution. 

In comparison to the observed velocity functions it is 
worth noting that the conversion between the line-of-sight 
stellar velocity dispersions and maximum circular velocity 
is somewhat uncertain (Stoehr et al. 2002; Zentner & Bul- 
lock 2003; Kazantzidis et al. 2003). Thus, at this point 
it is not worth trying to reproduce the observed function 
exactly. 

The median distances of luminous satellites to their re- 
spective hosts are 59ft. -1 kpc, 91/i _1 kpc, and 73/i _1 kpc 
for halos Gi, G2, and G3, respectively. This is close to the 
median values for the MW and M31 and lower than the 
median distance of the overall DM satellites (~ 120/i _1 
kpc, see § 5.4). The KS probability that the radial dis- 
tribution of these luminous satellites is drawn from the 
same distribution as that of the MW are 97%, 1%, and 
75% for the three hosts, respectively. Although there are 
apparent fluctuations due to the differences in the evolu- 
tionary histories of the three hosts, the luminous satellites 
in our model have a clear tendency to be more centrally 
concentrated than the overall DM satellite population. 

On the other hand, tidally-triggered bursts of star for- 
mation are not limited to the central parts of the host halo. 
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Fig. 9. — Morphological segregation of galaxies in the Local Group 
(lower panel) and in our model (upper panel). We divided the ob- 
served galaxies in two broad classes: Irr - all dwarf irregular galax- 
ies, and Sph - all spheroidal systems, including dSph, dSph/dlrr, 
and dEs. The model result include all subhalos that formed stars 
and the division into irregular and spheroidal systems was done us- 
ing the amount of heating experienced by each object, as explained 
in § 6.2. The model reproduces the observed morphological segre- 
gation: most spheroidal systems are located within 300ft" 1 kpc of 
the host, while irregular systems are found at a wide range of radii. 



In the isolated halo Gl, where the sample of satellites is 
not contaminated by the close proximity to another host, 
the tidal heating parameter / t id scales with the present 
distance as Jtid oc r°, a — —3.7 ± 0.2, for r < 1 h^ 1 Mpc, 
where the quoted error is the standard deviation of the 
whole sample, not the error of the mean. This is con- 
sistent with the expected slope a = —4 (c.f. eq. [A6], 
ignoring the adiabatic correction). However, if we limit 
the sample to the halos of interest, i.e. only those capable 
of forming stars (with the maximum V, n > V4 = 16.7 km 
s -1 , the virial velocity corresponding to T v i r = 10 4 K), 
then the slope is significantly shallower: a = —1.8 ± 0.2. 
Furthermore, if we consider only the largest tidal param- 
eters that might lead to starbursts (Itid > 10 3 Gyr~ 2 ), 
then the distribution is almost independent of distance: 
a = —0.3 ±0.1. Thus, the current location of the satellite 
in the host galaxy gives very little indication whether it 
had tidal starbursts in the past. 

We find a large variety of star formation histories for 
the luminous satellites. Most systems have a single ini- 
tial burst lasting up to 2 Gyrs. For some this is the only 
starforming activity, while other have a constant SFR at 
1 Mq yr" 1 up to z — 0.7 or bursty star formation continu- 
ing until z = 0.2. There are also objects that have only a 
single tidally-triggered burst at z ~ 1. The tendency is for 
more massive satellites to have more extended star forma- 
tion. The median mass-weighted epoch of star formation 



11 



in halos with V m < 30 km s~ x is between i moc i = 1 — 4 Gyr 
cosmic time (corresponding redshifts z = 5 — 1.7), while 
the more massive halos have £ mcc i up to 7 Gyr (z = 0.7). 

Our simple model also predicts central stellar densities 
in a reasonable agreement with observations: roughly con- 
stant £* ~ 5 - 50 M pc~ 2 for M* < 10 9 M and rising 
with the stellar mass as £* <~ M*/(10 7 M Q ) M Q pc~ 2 for 
systems with M* > 10 9 M©. The satellites located within 
100 kpc of their host galaxy have typically higher central 
densities (> 50 M© pc~ 2 ) than the more distant satellites. 

The results listed above are valid for all dwarf satellite 
galaxies regardless of their evolutionary history. In addi- 
tion, as discussed in the previous section, our model tracks 
the tidal heating of stars formed in each halo. We can 
therefore attempt a crude morphological classification of 
galaxies based on the amount of heating they experienced. 
This is motivated by the observations that dSph galaxies 
have low values of the ratio of rotation velocity to the ran- 
dom velocity dispersion, v ro t/a < 1. The galaxies of tran- 
sition type dlrr/dSph have v TOt /cr < 2 (see, e.g., Grebel 
et al. 2003, and references therein). Theoretical models 
of Mayer et al. (2001b, a) also indicate that the tidally- 
heated dSph-like remnants of low-surface brightness spiral 
galaxies have small v TOt jo. 

We use the circular velocity at the radius enclosing all 
of the stellar mass, v° ut , as a proxy for v TOt . The rota- 
tion velocity will, in general, be smaller than the circular 
velocity because some of the kinetic energy is in the form 
of the random motions. Also, we account only for direct 
tidal heating and do not take into account tidally-induced 
heating via bar and bending instabilities. The exact value 
of Wrot/c for our galaxies is thus somewhat uncertain as 
our a may be regarded as lower limit. We experimented 
with several values for the classification threshold in the 
range 1 < v° ut /<r < 3, but the main trends are not sen- 
sitive to the specific choice in this interval. We chose the 
value of u° ut /(j = 3 for the classification shown in Figure 9. 
For the observed galaxies, we combined dwarf spheroidal, 
transition type, and dwarf elliptical galaxies in one broad 
class of spheroidal systems, using Table 1 of Grebel et al. 
(2003). The figure shows that our model is consistent with 
the observed trend of morphological segregation. Most 
spheroidal systems are located within 300/i _1 kpc of the 
host, while irregular systems are found at a wide range of 
radii. The two model spheroidal systems at ~ l/i -1 Mpc 
have been part of a small group of galaxies and were tidally 
heated within this group before it was accreted by the host. 

7. DISCUSSION 

In the previous sections we presented the results of the 
dynamical evolution of galactic satellites in self-consistent 
cosmological simulations. One of the main findings is 
that the internal structure of the satellites responds to 
the changes of mass in a remarkably regular way. Namely, 
both during the periods of mass growth and tidal mass 
loss, the maximum circular velocity of a halo changes as 
U m cx M x t a . The slope a is ~ 4 — 5 when the mass grows 
and a w 3 when the mass decreases due to tidal stripping. 
The latter result was also obtained by Hayashi et al. (2003) 
in their non-cosmological simulations of satellite evolution. 

The overall evolution of subhalo population is such that 
their M — V m relation is similar to that of isolated halos 



(with a w 3.3). The fact that isolated halos and subhalos 
have similar mass-circular velocity relations may hint at 
why the fundamental plane of galaxies in clusters and the 
field arc similar (Dressier et al. 1987; Djorgovski & Davis 
1987; Mobasher et al. 1999; Bernardi et al. 2003) and why 
the scatter in the Tully-Fisher relation is so small (Kan- 
nappan et al. 2002). 

The fact that the circular velocity decreases with de- 
creasing mass means that the systems experiencing dra- 
matic mass loss will experience a significant change in 
circular velocity. We find that about 10% of the subha- 
los with masses < 10 8 - 10 9 M© or V m < 30 kms" 1 at 
z = have considerably larger masses and circular veloc- 
ities at earlier epochs. This may explain how such appar- 
ently small objects like Ursa Minor and Draco could have 
formed stars, given that the gas accretion is expected to 
be strongly suppressed by the UV background (e.g., Thoul 
& Weinberg 1996; Gnedin 2000). In our model, these sys- 
tems were once sufficiently massive (V m > 30 kms -1 ) to 
accrete gas and form stars but the accretion was halted 
when they started to experience tidal mass loss. 

After the accretion of new gas stops, these systems may 
continue to form stars in bursts as they are tidally stirred 
(e.g., Mayer et al. 2001b). Interestingly, we find that 
the strongest tidal interaction may occur even before halo 
is accreted by the host. Some satellites experience the 
strongest tidal force from multiple halos at early epochs 
in major mergers during the assembly of their host (see 
Fig. 3). Such mergers are frequent at early epochs, and we 
find that in general all satellites forming stars experience 
multiple bursts in the first 2 — 3 Gyrs of their evolution. 
We present a simple model for star formation in dwarf ha- 
los and apply it to the evolutionary tracks extracted from 
the simulations. As shown in Figure 7, the model is suc- 
cessful in reproducing the abundance of luminous satellites 
around M31 and the Milky Way. 

The spatial distribution of dwarfs around the Milky Way 
offers another independent challenge to any model of satel- 
lite evolution. Figure 8 shows that our model reproduces 
the observed distribution reasonably well. The distribu- 
tion of luminous satellites is more compact than the overall 
population of subhalos because stars form only in objects 
that were sufficiently massive at high redshifts. Due to the 
strong mass- and redshift-dependence of spatial bias, such 
objects are considerably more clustered around the host 
than smaller halos that form at a wide range of redshifts. 
Correspondingly, we find that luminous objects were ac- 
creted by the host systematically earlier (by Az ~ 0.5 — 1) 
than smaller mass dark subhalos. 

One of the remarkable features of our model is that the 
results are not sensitive to the details of reionization his- 
tory of the Universe. For example, all of the presented 
results are nearly intact if we change the assumed redshift 
of reionization from the fiducial value of z T = 7 to z r = 15 
(see § B). The physical reason behind this insensitivity 
to reionization is the inefficiency of gas cooling and star 
formation in small mass (T v i r < 10 4 K) systems. This 
is because gas in such systems cannot cool via hydrogen 
line emission and must rely on the inefficient H2 cooling. 
Such redshift-independent suppression of gas cooling is ob- 
served in cosmological simulations of Chiu et al. (2001) and 
Kravtsov & Gnedin (2003). The important implication is 
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that properties of the population of galactic satellites are 
determined by the physics of galaxy formation rather than 
by the UV background and reionization. 

Our results can qualitatively explain the morphologi- 
cal segregation of the Local Group galaxies (e.g., Grebel 
2000). As shown in Figure 9, a simple division of model 
galaxies into irregular and spheroidal based on the amount 
of tidal heating they experienced during their evolution 
reproduces the main observed trend. Most spheroidal sys- 
tems are located within 300ft. -1 kpc, while irregular galax- 
ies are found almost uniformly at all distances. Therefore, 
our results support the scenario that spheroidal systems 
form via strong tidal heating (Mayer et al. 2001b, a). Note, 
however, that tidal heating is not restricted to the host. 
It can occur early on, before the host is assembled, within 
merging subgroups. 

Interestingly, this explains a puzzling presence of the 
Cetus and Tucana dSph galaxies at the outskirts of the 
Local Group some 700 kpc and 1000 kpc from the nearest 
massive spiral (MW or M31). We also find ^1 — 2 galaxies 
with significant heating (the ratio of the rotational veloc- 
ity to the velocity dispersion of v TOt ja < 1) at distances 
<~ 1000 kpc from their hosts. The tidal heating of these 
systems occurred in small groups that are being accreted 
by the host at the current epoch (see also Gnedin 2003b, 
for a similar effect in clusters of galaxies). As the tidal 
force unbinds satellites from such accreting groups, iso- 
lated dSph galaxies may be found at large distances from 
the primary host. 

Also, early tidal interaction, experienced for example by 
the system shown in the left column of Fig. 3, and sub- 
sequent interaction with other subhalos may lead to the 
increase of orbital energy and apocenter distance. This 
scenario would also explain presence of dSphs at large dis- 
tances from the primary. The main point in both scenar- 
ios is that primary is not the only source of tides and the 
present-day environment is not necessarily indicative of a 
dwarf galaxy's past. 

One of the most interesting candidates for the "miss- 
ing" dark halos is the population of compact high-velocity 
clouds (CHVCs) of neutral hydrogen (HI, e.g., Blitz et al. 
1999; Braun & Burton 1999). This idea has recently been 
boosted by the detection of concentration of CHVCs near 
M31 (Thilker et al. 2003). It is thus interesting to con- 
sider the amount of gas associated with the subhalos that 
remain dark in our model. The cumulative gas mass func- 
tion associated with dark halos is remarkably consistent 
for all three host halos: N(> M g ) « 20(M g /10 7 M )-°- 7 
for 10 6 < M g < 10 8 M within 200ft" 1 kpc. Most of 
the gas mass is thus in most massive subhalos. The total 
mass of gas associated with such halos within 200ft~ 1 kpc 
is Mg 0t « 2 x 10 9 M Q , the number similar for all three 
hosts. If we assume that on the average about 10% of gas 
is neutral (Maloney & Putman 2003; Thilker et al. 2003), 
the total mass in neutral hydrogen is Mhi - 2x 10 8 M Q . 

The predicted number of dark clouds with M g > 10 6 M Q 
is ~ 50 — 100. A fraction of the observed CHVCs can thus 
be associated with the small-mass DM halos. Within cen- 
tral 50 kpc, however, the number of halos with such gas 
masses is only ^2 — 5. We cannot therefore explain 25 
CHVCs observed by Thilker et al. (2003) within this radius 
around M31. It is possible that simulations underpredict 



the number of small-mass halos due to overmerging. To 
check this will require higher- resolution simulations. On 
the other hand, we did not take into account processes 
such as ram pressure stripping, which would further re- 
duce the number of halos with gas. Another possibility is 
that most of the observed M31 CHVCs are gas clouds in 
tidal streams, such as the Magellanic Stream, and are not 
associated with distinct dark matter halos (Putman et al. 
2003). 

8. COMPARISON WITH PREVIOUS WORK 

Possible astrophysical solutions 6 to the "missing satel- 
lite problem" have been considered in the last several years. 
Here we discuss the main differences of our model and the 
models proposed in previous studies. 

Bullock et al. (2000), Somerville (2002), and Benson 
et al. (2002) discussed the formation and evolution of dwarf 
galactic satellites using semi-analytic models of different 
degrees of sophistication. The conclusion reached by all 
these studies is that the extragalactic UV background can 
greatly suppress the gas accretion and star formation in 
the majority of low-mass (V m < 30 kms -1 ) halos. A small 
fraction of the dwarf halos that harbors stellar systems 
was assumed to have formed (i.e., assembled significant 
fraction of their mass) before reionization, when the level 
of UV radiation was low. This is because in all of these 
studies the maximum circular velocity of subhalos was as- 
sumed to be constant as the mass is tidally stripped. There 
was thus a simple one-to-one mapping between the circu- 
lar velocity observed at z = and at the time of accretion. 
Our results show that this assumption is incorrect (see also 
Hayashi et al. 2003; Kazantzidis et al. 2003). Another key 
difference is that tidal mass loss in our model can occur 
before a halo is accreted by the host, as a result of inter- 
actions with other halos. These effects are not accounted 
for in any of the semi-analytic models. 

The implicit assumption in the above models is that the 
small systems would be able to retain the accreted gas and 
form stars after reionization. This assumption was justi- 
fied at the time, as the first calculations of photoevapora- 
tion of gas indicated that halos with V m > 10 km s _1 might 
retain their gas (Barkana & Loeb 1999). More recent cal- 
culations, however, show that the gas could be gradually 
removed from halos of up to V max ~ 30 kms -1 (Shaviv & 
Dckcl 2003). In light of this result, the previous models 
would not be able to explain the formation and properties 
of luminous dwarfs, as the star formation in small halos 
would be suppressed after reionization. It would thus be 
difficult to explain the more extended star formation his- 
tories derived for many dSph galaxies in the Local Group 
(Grebel 2000). 

In our model, the small-mass dwarfs are identified with 
the halos that were relatively massive at high redshift and 
could retain the gas and form stars after reionization. The 
star formation histories of dwarfs are thus more extended, 
in better accord with observations. As noted in the previ- 
ous section, our model is also insensitive to the epoch of 
reionization and can accommodate early reionization sug- 

6 The solutions that invoke astrophysics of galaxy formation within 
the standard CDM framework rather than modifications of the prop- 
erties of dark matter particles and/or the shape of the initial power 
spectrum. 
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gested by results of the WMAP satellite (Spcrgel et al. 
2003). 

Our model and all of the models discussed above are 
qualitatively different from the proposal of Stoehr et al. 
(2002, 2003). These authors argued that the maximum cir- 
cular velocity of the Local Group dwarfs may be systemat- 
ically underestimated because it is derived from the stellar 
velocity dispersion within radii considerably smaller than 
r max , the radius at which the maximum halo velocity, V m 
is reached (see, however, Kazantzidis et al. 2003). Stoehr 
et al. (2002, see also Hayashi et al. (2003)) suggested that 
the luminous dwarfs may be harbored by the most massive 
satellites of the DM halos. This has an important physical 
implication: if the dwarfs indeed occupy twelve or so most 
massive halos, then there exists a certain mass scale below 
which galaxy formation is completely suppressed. If, on 
the other hand, the dwarf galaxies occupy satellites with 
a variety of masses (<~ 10 7 — 10 10 M Q ), one has to explain 
why some fraction of small halos managed to light up the 
stars, while most others did not. 

If the idea of Stoehr et al. (2002) is correct, our re- 
sults indicate that circular velocities of dwarf spheroidal 
halos should have been even larger (by a factor of two or 
more) than the values inferred from the current observa- 
tions. This could make halos of some galaxies uncomfort- 
ably massive. For example, Stoehr et al. (2002) derive 
the maximum circular velocity for the Draco in the range 
~ 35 — 55 kms -1 . This implies the pre-accretion values of 
Knax 70 kms -1 and the pre-accretion mass comparable 
to those of M32, NGC 205, and M33. The fact that lumi- 
nosity of Draco is almost four orders of magnitude lower 
than luminosities of these galaxies would present a major 
puzzle. 

In addition, the radial distribution of the most massive 
satellites should be consistent with the observed radial dis- 
tribution of the MW satellites. We find that in our sim- 
ulations the radial distribution of subhalos with largest 
V m is between that of the luminous satellites and all DM 
satellites shown in Figure 8. In a study of a larger sam- 
ple of cluster halos, De Lucia et al. (2003) find that the 
radial distribution of the most massive halos is even more 
extended than that of the smaller mass objects. A similar 
point was made recently by Taylor et al. (2003), who used 
semi-analytic models for subhalo population to show that 
the radial distribution of the most massive halos is more 
extended than that of the MW satellites at > 3cr level. A 
caveat to this argument is that the sample of Milky Way 
satellites may be incomplete at large distances and more 
faint dwarf galaxies will be discovered in the future (Will- 
man et al. 2004). 

9. CONCLUSIONS 

We presented a study of the dynamical evolution of 
galactic satellites using self-consistent high-resolution cos- 
mological simulations of three MW-sized halos. Our main 
results and conclusions are as follows. 

• We find that ss 10% of the substructure halos that 
have masses of < 10 8 — 10 9 M at the present epoch, 
had considerably higher masses and circular veloc- 
ities when they formed at z > 2. After the initial 
period of mass accretion, while these objects evolve 



in isolation, they suffer dramatic mass loss due to 
tidal stripping by actively merging massive neigh- 
boring halos. Strong tidal interactions can occur 
even before the dwarfs are accreted by their pri- 
mary host halos. 

• The decrease in mass due to tidal stripping is ac- 
companied by the decrease in the maximum circu- 
lar velocity, such that the objects evolve along a 
M - V* relation with a w 3 - 4. 

• These results indicate that some of the systems that 
have small masses and circular velocities at z = 
could have had masses comparable to those of the 
SMC and LMC in the past. This can explain how 
the smallest dwarf spheroidal galaxies observed in 
the Local Group were able to build up sizable stellar 
masses in such shallow potential wells. 

• We present a simple galaxy formation model based 
on the evolutionary tracks extracted from the sim- 
ulations. The novel features of the model are the 
starburst mode of star formation after the strong 
peaks of the tidal force and accounting for the inef- 
ficient dissipation of gas in halos with T vlr < 10 4 K. 
The model can successfully reproduce the circular 
velocity function, radial distribution, morphological 
segregation of the observed Milky Way satellites, 
and the basic properties of galactic dwarfs such as 
stellar masses and densities. 
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APPENDIX 
CALCULATION OF THE TIDAL FORCE 

In our analysis we use the external tidal force experi- 
enced by each satellite halo to estimate the strength of 
tidal interaction. We calculate the force both directly from 
the gravitational potential computed in the simulation and 
using an analytical approximation for the neighbor halos. 

To compute the tidal force numerically from the local 
potential we estimate its second spatial derivative at 
the center-of-mass of the satellite: 

F «=-{d^X r ^ Fa ^ (A1) 

where r is the radius- vector in the satellite reference frame 
and R is the radius- vector in the perturbcr reference frame. 
The potential $ is calculated on the original refinement 
grid using the ART gravity solver. In the calculation of 
the potential, we subtract the self contribution of the halo 
and consider only the external tidal potential. 

In a study of galaxy interactions in clusters of galaxies, 
Gnedin (2003b) used the Savitzky-Golay smoothing filter 
to interpolate the potential on a plane and calculate its 
derivatives from a smooth polynomial function. We cm- 
ploy a similar scheme but with the adaptive 4-th order 
interpolating polynomials in each of the three orthogonal 
planes around the satellite center of mass: 

4 

P 4 (x,y) = c kix k y l (A2) 

fc,Z=0 

and the same for the xz and yz planes. The 4-th order 
expansion ensures a smooth second derivative of the po- 
tential. In each of the planes we extract a n x n subgrid 
centered on the original grid point, nearest to the satellite 
center. In order to obtain a uniform accuracy of the tidal 
force for satellites of different sizes, we choose the size of 
the subgrid cells to be closest to 1/4 of the satellite's tidal 
radius. The coefficients c^i are calculated by minimizing 
X 2 deviation 

n 

x 2 = ^ x ^vi) - ^yt)f ( A3 ) 

using the CERN Program Library routine MINUIT 7 . We 
have experimented with n = 16, 32, and 64 and found that 
n = 64 provides the most accurate derivatives, as tested 
on the analytical NFW models. The tidal tensor compo- 
nents F a p are then calculated by analytical differentiation 
of equation (A2). 

7 http : //wwwasdoc . web . cern . ch/wwwasdoc/minuit/ 
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We compare the real tidal force due to the overall mass 
distribution in the simulation with the contributions of all 
neighboring halos, including the host halo. We model the 
halos with an NFW density profile and take their mass 
M v ; r and virial radius r v j r directly from the halo cata- 
logs generated by the halo finder (see § 3). We determine 
the scale radius of the NFW model for the satellite ha- 
los from the position of the maximum circular velocity, 
r s = f max/2. 16. For the host halo, we use the parametriza- 
tion c n f w = r v i r /r s = 16a 3 / 2 , which is a best fit to the 
density profile of the analyzed host halos. The analytical 
tidal force in the reference frame of the satellite is then 
readily calculated using eq. (5) of Gnedin et al. (1999): 
GM(R) 



F(r) = 



R 3 



[(3-M)(n-r)-r] 



(A4) 



where r is the radius-vector within the satellite, R is the 
distance to the perturber, n = R/i?, (i = rflnM /d\nR, 
and M(R) is the enclosed mass of the NFW model: 
In (1 + R/r s ) - 1 + (1 + R/rs)- 1 



M{R) = M v 



(A5) 



ln(l + C„f w ) - 1 + (1 + c n f w )- ] 
Figure 3 shows that the approximate tidal force calculated 
in this manner is quite accurate, especially near the max- 
imum of the tidal force. 

Although the tidal force along the satellite trajectory 
varies rapidly with time, most of the tidal heating of stars 
and dark matter particles occurs near the strong peaks of 
the tidal force. Each of these tidal peaks can be considered 
as an independent tidal shock (Gnedin & Ostriker 1999; 
Gnedin 2003b) . The amount of tidal heating, such as the 
increase of the velocity dispersion, is proportional to the 
integral over the peak of tidal force: 

-3/2 

72^ I , (A6) 




Itid(t n ) 



where the sum extends over all components of the tidal 
tensor, a,(3 = {x,y,z}. The last factor is the correction 
for the conservation of adiabatic invariants of stellar or- 
bits during the tidal shock (c.f., Gnedin & Ostriker 1999). 
Here r„ is the effective duration of peak n at time t n , 
and tdyn is the dynamical time of the satellite. We take 
idyn = 27rr 1 /2/t>rot , where 7^/2 is the half-mass radius of 
the stellar disk and v mt is the circular velocity of the appro- 
priate NFW model at r 1 / 2 - The cumulative tidal heating 
parameter is the sum over all tidal peaks: 

/tid = ^itid(U- (A7) 

n 

This parameter determines the increase of the velocity dis- 
persion of stars (eq. [6]) in our model of dwarf galaxy 
formation (§6). 

FILTERING MASS SCALE 

We estimate the suppression of gas accretion due to 
the extragalactic UV background using the filtering mass 
scale, derived by Gnedin (2000). He defined Mp as the 
mass of the halo which would lose half of the baryons, 
compared to the universal baryon fraction. This filtering 
mass relates to the Jeans mass of the intergalactic gas in- 
tegrated over the cosmic history (eq. [6] in Gnedin 2000): 
M F (a)=M J0 /(a) 3 / 2 , (Bl) 
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Fig. B10. — Top panel: filtering virial circular velocity (V c ) and 
virial temperature (axes are chosen such that the two solid lines co- 
incide), and maximum circular velocity (Vm, dashed line). Bottom 
panel: filtering mass. Data points with error-bars show the sim- 
ulation results of Gnedin (2000). Solid lines are for the standard 
epoch of rcionization: z r = 7, z = 8. For comparison, the dotted 
line in the bottom panel shows the case of earlier reionization with 
z r = 10, z = 11. 



where M J0 = 2.5 x lO 11 /*- 1 ^,, 1/2 /i" 3/2 M Q , /x w 0.59 is 
the mean molecular weight of the fully ionized gas, and 
the integration extends over the expansion factor, a. The 
temperature of the cosmic gas T4 is expressed in units of 
10 4 K for convenience. 

Here we propose an analytical fit to the results of Gnedin 
(2000), assuming a simple dependence of the temperature 
on the expansion factor: T^ia) = (a/a ) a for a < a , 
T^{a) = 1 for a Q < a < a r , and T^a) = (a/a r ) _1 for 
a > a r . 

These three distinct stages can be clearly seen on Fig. 
1 of Gnedin (2000). They correspond, to the 1) epoch 
before the first HII regions form, z > z D , 2) the epoch of 
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the overlap of multiple HII regions, z r < z < z a , and 3) 
the epoch of complete reionization, z < z r . In the first 
stage, before redshift z Q = l/a — 1 w 8, the temperature 
is rising as the newly-formed stars ionize their neighboring 
regions. The parameter a controls the rate of growth of 
the extragalactic UV flux; we find a = 6 to be the best fit. 
During the overlap stage, between redshifts z Q and z r = 
1 /a r — 1 w 7, the temperature is kept constant at roughly 
10 4 K as the cosmic HII regions overlap. After the universe 
is fully ionized, at redshifts below z r , the temperature falls 
adiabatically with the cosmic expansion. 

With these analytical expressions for T^a), we integrate 
equation (Bl) analytically: 
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The virial circular velocity of the halo is 
V c 3 = GMiJ(z)(A vir /2) 1 /2 j where = Jf(0)[ft (l + 

z) 3 + ilA] 1//2 is the Hubble constant, and A v ; r is the virial 
overdensity with respect to the critical density, parametrized 
by Bryan & Norman (1998) as A vir (z) = 18tt 2 + 82.x - 
39x 2 , x = Q(z) — 1. The virial temperature is T v ; r = 
36(V c /km s" 1 ) 2 K. 

Our analytical fit is convenient for accurate modeling of 
the photoheating effect in semi-analytical models of galaxy 
formation. Its versatile form, with two parameters z and 
z r , allows a simple recalculation of the filtering mass for 
a different redshift of reionization than was assumed in 
the simulation of Gnedin (2000). It can also be easily 
adapted to describe two epochs of reionization, or the early 
extended reionization suggested by WMAP (Spergel et al. 
2003). For illustration, we show on the lower panel of 
Figure B10 the filtering mass as a function of redshift for 
two choices of the reionization redshift. The top panel 
shows the filtering circular velocity and the corresponding 
values of the maximum velocity V m and virial temperature 
T ■ 

± vir ■ 

Gnedin (2000) provided the following expression for the 
amount of cold gas left in the halo of mass M: 

^■^ [i + oAvmP - (B3) 

where /b ~ 0.14 is the universal baryon fraction. In § 6 
we use this fraction of cold gas to model the star forma- 
tion history of the satellite galaxies. To account for the 
inefficiency of atomic gas cooling at T < 10 4 K, we apply 
equation (B3) substituting for Mp the maximum of M-p(z) 
or M4, the halo mass corresponding to T v ; r = 10 4 K. 



